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ALTERNATING DIRECTION IMPLICIT METHODS FOR PARABOLIC 
EQUATIONS WITH A MIXED DERIVATIVE 

RICHARD M. BEAM + and R. F. WARMING 1 " 

Abstract . Alternating direction implicit (ADI) schemes for two- 
dimensional parabolic equations with a mixed derivative are constructed by 
using the class of all A Q -s table linear two-step methods in conjunction with 
the method of approximate factorization. The mixed derivative is treated with 
an explicit two-step method which is compatible with an implicit A c -stable 
method. The parameter space for which the resulting ADI schemes are second- 
order accurate and unconditionally stable is determined. Some numerical 
examples are given. 


CONTENTS 

1. Introduction 2 

2. Preliminaries 5 

2.1. Linear multistep methods 5 

2.2. One-leg (multistep) methods .... 8 

2.3. Combined linear multistep methods 9 

3. Unconditionally stable schemes 10 

4. An ADI scheme: the p(E) or A formulation 15 

5. A general two-step ADI scheme 19 

5.1. General formulation 19 

5.2. Special cases of general formulation 21 


4 . 

'Computational Fluid Dynamics Branch, Ames Research Center, NASA, Moffett 


Field, California 94035. 



5.3 Algorithm selection 24 

5.4 General formulation with no mixed derivative 26 

6. Time-dependent coefficients 27 

7. Numerical examples 29 

8. Concluding remarks 34 

Appendix A: Stability analysis of combined LMMs 36 

Appendix B: Stability analysis for two-step ADI schemes ... 38 

References 43 


1. Introduction. When an alternating direction implicit (ADI) method is 
applied to a parabolic equation, for example, 

(1.1a) du(yi^y,t) _ Lu ( x>yjt ) # 


where 

g2 g2 g2 

(1.1b) L = a(x,y,t) + b(x,y,t) + c(x,y,t) j-j , 

(1.1c, d) a,c > 0 , b 2 < 4ac , 

it reduces the computational problem to a sequence of one-dimensional (matrix) 
inversion problems. If the mixed derivative 3 2 /3x3y of the operator L 
were absent (b = 0) , an ADI method would reduce the operator (1 - L) to the 
product of two one-dimensional spatial operators. In the method of Douglas 
and Gunn [9], the mixed derivative is kept implicit and their scheme requires 
four inversions; that is, the operator (1 - L) is reduced to four one- 
dimensional operators. The two additional inversions are the direct result of 
keeping the mixed derivative implicit. Simpler and more efficient schemes 
can be obtained if the mixed derivative is evaluated explicity. Perhaps 
not so obvious is the stability of these simpler schemes. In fact, it is 
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rather surprising that one can develop unconditionally stable algorithms for 
(1.1) by computing the mixed derivative explicitly and the derivatives 
8 2 /3x 2 and 3 2 /3y 2 implicitly. McKee and Mitchell [15] surveyed two-level, 
first-order accurate (in time) schemes and devised a new first-order accurate, 
unconditionally stable scheme for (1.1). Iyengar and Jain [12] generalized 
the method of McKee and Mitchell and presented a three-level, second-order 
accurate scheme for (1.1); however, the implementation of the scheme is com- 
plicated by the explicit computation of fourth differences although the origi- 
nal partial differential equation (1.1) contains only second derivatives. 

The present authors constructed a second-order accurate ADI algorithm for 
the compressible Navier-Stokes equations [1]. When the algorithm is applied 
to the model equation (1.1), it leads to a simple three-level unconditionally 
stable scheme which is easy to implement. In a recent paper [20] we combined 
A-stable linear multistep methods (LMMs) and approximate factorization to 
construct a large class of multilevel unconditionally stable ADI schemes for 
a model partial differential equation with both convective (hyperbolic) and 
diffusive (parabolic) terms; that is, 

T _ 3 2 , . 3 2 , _3£_ _3_ _ _3_ 

L 3x 2 b 3x3y + C 3y 2 C 1 3x C 2 3y ’ 

where a,b,c satisfy (1.1c, d) and c 1# c 2 are real constants. The general 
formulation of [20] is second-order accurate if b = 0 but first-order accu- 
rate in time if the mixed derivative is included. In both [1] and [20], the 
mixed derivative is treated explicitly. The purpose of the present paper is 
to modify and combine the algorithms of [1] and [20] to obtain a general, 
second-order accurate, unconditionally stable algorithm for the model parabolic 
equation (1.1). In a companion paper [21] we apply the method to derive a 
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noniterative ADI algorithm for a hyperbolic-parabolic system of nonlinear 
equations with mixed derivatives. 

The development and analysis of numerical methods for ordinary differen- 
tial equations (ODEs) is more advanced than that for partial differential 
equations (PDEs). Therefore, it seems plausible to capitalize on this fact by 
making use of known results from the theory of difference methods for ODEs to 
construct methods for PDEs. For example, the time differencing schemes used 
to construct implicit methods for PDEs are invariably LMMs although this fact 
is seldom noted. Since a great deal is known about the properties of LMMs 
(see, e.g. , [6,10]), one can use this information to advantage when attempting 
to construct schemes for PDEs. With these observations in mind we use Sec. 2 
as a review of the theory and notation for linear multistep methods including 
A-stability, A 0 -stability, and one-leg methods. In addition, we introduce the 
notion of an LMM combining two different LMMs — one implicit and the other 
explicit. In Sec. 3 we investigate the stability of an implicit method for 
(1.1) obtained by using an A Q -stable LMM as the time differencing method. We 
then modify the scheme by treating the mixed derivative explicitly and rein- 
vestigate the stability properties. The method of approximate factorization 
is applied in Sec. 4 to obtain an unconditionally stable, second-order accurate, 
multistep ADI scheme in the p(E) formulation. In Sec. 5 we construct a gen- 
eral approximate factorization method for all second-order, two-step schemes 
and discuss the parameter space for unconditional stability. Details of the 
stability analyses are given in the Appendices A and B. In Sec. 6 a simple 
modification is given for the case of time- dependent coefficients a,b,c. 
Numerical examples are given in Sec. 7 and some concluding remarks in Sec. 8. 
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2. Preliminaries. In this section we briefly review the theory of 
linear multistep methods (LMMs) and the related one-leg (multistep) methods. 
In addition, we introduce combined LMMs. 

2.1. Linear multistep methods. A linear k-step method for the first- 
order ordinary differential equation 

du 


— = f(u,t) , u(0) = u Q 


( 2 . 1 ) 

is defined by the difference equation 
( 2 . 2 ) 

where p and 0 are the generating polynomials, 


p(E)u n = Ato (E) f n , 


(2 . 3a, b) 


p(o = 22 Y j , o(c) - 22 B j ?:i » 

j=0 

and E is the shift operator, that is. 


j = 0 


(2.4) 


„ n n+1 
Eu = u 


In (2.2), u n is the numerical solution at the point t = t n = nAt, At is the 
time step, and f n = f(u n ,t n ). The method is explicit if = 0 and implicit 
otherwise. Consistency and normalization are expressed by the relations: 

k k k 

(2.5a,b,c) 22 °j = 0 » 2D ja j = S 6 j = 1 • 

j = 0 j = 0 j = 0 

As an example of an LMM, the most general consistent two-step method 
(i.e., k = 2 in (2.3)) can be written as 

(2.6) (1 + C)u n+2 - (1 + 2£)u n+1 + £u n = At[0f n+2 + (1 - 0 + <{>)f n+1 - <(>f n ] 

where (0,^,<}>) are arbitrary real numbers. The operators p(E) and o(E) are 
(2.7a) p(E) = (1 + C)E 2 - (1 + 25)E + E, , 

(2.7b) a(E) = 0E 2 + (1 - 0 + <f>)E - 4> . 
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For the class of all two-step methods that are at least second-order accurate, 
the parameters (0,5,<}>) are related by 


(2.8) <j> = K ~ 6 + \ , 

and consequently, cr(E) can be rewritten in terms of the two parameters (0,5) as 

(2.9) o(E) = 0E 2 + (5 - 26 + |)e - (g - 0 + ) . 

Some well-known implicit second-order methods and their corresponding values 
(0,5,<j>) are listed in Table 1. Linear one-step methods are a subclass of (2.6) 
obtained by setting 5 = <J> = 0: 

(2.10) u n+1 - u n = At [0f n+1 + (1 - 0)f n ] , 

where we have shifted the time index down by one. The trapezoidal formula 
(0 = 1/2 in (2.10)) 

(2.11) u - u = ~ 2 ~ (f +f) 


is the only second-order accurate one-step method. When the trapezoidal for- 
mula is applied to a parabolic equation, the resulting algorithm is usually 
called Crank-Nicolson. 

The linear stability of an LMM is analyzed by applying it to the linear 
test equation 


( 2 . 12 ) 


du 

dt 


= Xu 


where X is a complex constant. The stability is determined by the location 
of the roots of the characteristic equation, 

(2.13) p(£) - XAta(c) = 0 , 

relative to the unit circle in the complex plane. The stability region of an 
LMM consists of the set of all values of XAt for which the characteristic 
equation (2.13) satisfies the root condition; that is, its roots ^ all 
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satisfy < 1 and the roots of unit modulus are simple. An LMM is said 

to be A-stable if its stability region contains all of the left half of the 
complex AAt plane including the imaginary axis [4]. A simple test for 
A-stability can be formulated in terms of positive real functions [7]. By 
applying the test to the linear two-step method (2.6), one finds that the 
method is A-stable if and only if 


(2.14a) 

(2.14b) 

(2.14c) 


0 > * + j , 


£ > -— 

* - 2 ’ 


5 < 6 + <P ~ j 


An LMM is 
tive real 
Again, by 
(2.6) is 
(2.15a) 

(2.15b) 

(2.15c) 


said to be A Q -stable if the region of stability contains the nega- 
axis of the complex AAt plane, that is, the interval (-», 0], 
applying the theory of positive real functions, one finds that LMM 
A Q -stable if and only if 

6 > <f> + \ > 

£ > -± 

- 2 ’ 

0 < 9 + <J> . 


The details of the analysis for obtaining these inequalities are described in 
[21], For first-order accurate schemes, the inequality (2.15c) is less strin- 
gent than (2.14c). This is not surprising since A Q -stability is a weaker 
requirement than A-stability. However, for the class of all second-order 
methods, the parameters (0,£,<ji) are related by (2.8), and both sets of inequal- 
ities (2.14) and (2.15) reduce to 

(2. 16a,b) £ < 20 - 1 , K > . 


7 


The parameter space (0,5) for which the class of two-step, second-order methods 
is A-stable, happens to coincide with the parameter space for which the class 
is A 0 -stable and is shown by the shaded region of Fig. 1. The methods listed 
in Table 1 and indicated by the symbols in Fig. 1 are both A Q - and A-stable. 

Dahlquist has proved that the order of accuracy of an A-stable LMM 
cannot exceed two [4]. On the other hand, Cryer [3] has proved there exist 
A Q -stable LMMs of arbitrarily high order. Since the eigenvalue associated 
with the parabolic equation (1.1) is real and negative, the application of an 
A Q -stable LMM will yield an unconditionally stable scheme, that is, no stabil- 
ity restriction on the size of At. 

In this paper we restrict our attention to second-order accurate LMMs. 

This limitation is motivated by two practical considerations. First, conven- 
tional techniques for constructing alternating direction implicit (all 
A-stable and A Q -stable LMMs are implicit [4,3]) schemes generally impose a 
second-order- temporal accuracy limitation independent of the accuracy of the 
LMM chosen as the time differencing approximation. The reason for this will 
become apparent in Sec. 4. In principle, by altering conventional proce- 
dures, ADI schemes of temporal order greater than 2 can be constructed for 
PDEs where A Q -stable LMMs are appropriate; however, the unconditional stabil- 
ity of the higher order ADI schemes is an open question. The second practical 
consideration is computer storage. Even if the question of unconditional 
stability is answered in the affirmative, the scheme must be at least a three- 
step method since the two-step method (2.6) contains no A 0 -stable third-order 
subclass . 

2.2. One-leg (multistep) methods. A class of methods closely related to 
LMMs is the one-leg (multistep) methods. The one-leg (k-step) method [5] 
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corresponding to (2.2) is 


(2.17) p(E)u n = Atf(o(E)u n , a(E)t n ) 

where u n denotes the one-leg method solution. Formally, the one-leg method 
(2.17) can be obtained by shifting the operator a(E) inside the argument 
parenthesis of f n = f(u n ,t n ) in the linear multistep formula (2.2). As an 
example, the trapezoidal formula (2.11) is an LMM and the implicit midpoint 
rule, 


(2.18) 


^n+1 

u 


^n 

u 


An+1 
Atf (— ; 


+ u 


n 



is the corresponding one-leg method. Dahlquist [5] has proved the following 

r A II T 

theorem relating solutions of a one-leg method and an LMM: Let {u } be a 

vector sequence that satisfies the one-leg difference equation (2.17) and set 

u = a(E)u 

Then {u 11 } satisfies the LMM difference formula (2.2). 

It should be noted that LMMs and one-leg methods are equivalent when 
applied to the linear test equation (2.12) and, consequently, the results of 
linear stability analysis are the same for both. However, one-leg methods 
simplify nonlinear stability analysis [5] and, in addition, they are more 
reliable than LMMs when used with rapidly varying integration step sizes [17], 
As applied in this paper, the one-leg formulation makes it easy to construct 
an ADI scheme for the parabolic Equation (1.1) with time-dependent coeffi- 
cients . 


2.3. Combined linear multistep methods. If the function f (u) of the 
differential equation (2.1) is split into a sum, that is, 

(2.19) 4r = f(u,t) = f x (u,t) + f 2 (u,t) , 
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we can construct an integration formula by combining two different LMMs. For 
example, 

(2.20) p(E)u n = AtOj(E)f“ + Ato 2 (E)f 2 , 

where the subscripts on Oj and cr 2 indicate that the coefficients B j of the 
generating polynomial (2.3b) differ for the two functions f} and f 2 . The 
analogous split one-leg method (2.17) is 

(2.21) p(E)u U = Atf j (Oj (E)u 11 , aj(E)t n ) + Atf 2 (cr 2 (E)u 11 , c? 2 (E)t n ) . 


3. Unconditionally stable schemes. In this section we examine the sta- 
bility of an implicit method for the parabolic equation (1.1) where the time 
differencing approximation is an A 0 -stable LMM. Next we modify the scheme 
by treating the mixed derivative explicitly and determine the criteria for 
unconditional stability. Although not essential to the final goal of con- 
structing unconditionally stable ADI schemes with the mixed derivative com- 
puted explicitly, this intermediate analysis does isolate the stability con- 
straints due to the application of combined LMMs and those due to approximate 
factorization. Furthermore, the nonfactored scheme is formulated so that the 
ADI variant follows directly (Sec. 4). 

Numerical methods for solving the parabolic equation (1.1) can be obtained 
by a direct application of the LMM (2.2). Since our interest is in construct- 
ing unconditionally stable schemes, we assume that the LMM is A Q -stable, By 
comparing (1.1) and (2.1), we identify 


(3.1) 


f(u) = Lu = 





where L is a linear differential operator. For simplicity we assume that 
the coefficients a,b,c are independent of time. The case of time-dependent 
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coefficients is considered in Sec. 6. Insertion of (3.1) into (2.2) yields 


(3.2) p(E)„"- i t^ + b^ 7+ c^o(E)u'' . 

To analyze the stability of (3.2) we assume a solution for the PDE (1.1) 
(with constant coefficients) of the form 

(3.3) u(x,y, t) = v(t)e 1 ^ <lX+IC2y ^ 

where v(t) is the Fourier coefficient and k 2 are the Fourier variables 

(wave numbers). Substitution of (3.3) into (1.1) yields an ODE for v(t): 

(3.4a) = Xv , 

where 

(3.4b) X = -(aKj 2 + bK j k 2 + ck 2 2 ) . 

The quadratic form in the parenthesis of (3.4b) is positive definite if and 
only if the inequalities (1.1c, d) are satisfied. These constraints consti- 
tute the parabolicity condition of the PDE and ensure that the solution of 

(3.4) is damped with time. To complete the stability analysis of the PDE 
scheme (3.2), we need only consider the stability of the LMMs (2.2) applied to 

(3.4) with X < 0. Since we assumed that (2.2) is A c -stable, the PDE scheme 
(3.2) is unconditionally stable. In practice, the spatial derivatives in (3.2) 
are replaced by appropriate difference quotients; however, as shown at the end 
of Appendix B, central spatial discretization does not alter the unconditional 
stability criteria obtained by assuming spatially continuous solutions. 

For the one-step methods (2.10), the scheme (3.2) reduces to 

/-o o n+l n / 9 2 , , 3 2 , 3 2 \ r „ n+1 , ,, n, 

(3.5) u - u = At [a — r + b - + c — yj[0u + (1 - 6)u ] . 

y 3x 2 dx °y 3y v 

With central spatial difference approximations, (3.5) is identical to a 
scheme suggested by Lax and Richtmyer [13]. Their paper appeared about the 
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same time the original ADI methods were proposed by Peaceman and Rachford [18] 
and Douglas [8]. The first ADI methods [8, 18] did not include a mixed deriv- 
ative term and its presence precludes the construction of an efficient ADI 
method. A simple way of circumventing this difficulty is to treat the mixed 
derivative explicitly. One might expect, however, that this would have an 
adverse effect on the unconditional stability of the algorithm. This stability 
question is considered in the remainder of this section. 

Consider the combined LMM (2.20) where Oj and a 2 are defined as follows. 
Let (2.2) represent a second-order A Q -stable LMM and define 
(3. 6a, b) o 1 (E) = a(E) , a 2 (E) = a e (E) , 

where a e (E) is a second-order explicit LMM (i.e., 8=0 in the generating 

1C 

polynomial (2.3b)) with the same generating polynomial p(?) as for the 
A Q -stable LMM. With these definitions, (2.20) becomes 

(3.7) p(E)u n = Ata(E)fJ + Ata e (E)fJ 

Henceforth, in reference to a combined implicit-explicit method such as (3.7), 
we refer to the LMM that defines p(E) and a(E) as the generating LMM. The 
linear stability properties of (3.7) for second-order two-step methods are 
examined in Appendix A. 

For didactic purposes in this section and practical reasons in the fol- 
lowing section, we rewrite the LMM (3.7) as 

(3.8) p(E)u n - coAtp (E)f j 1 = At [a (E) - wp(E)]f“ + Ata e (E)f“ , 
where 

(3.9) OJ = . 

The parameter u> is defined so that the operator a(E) - top (E) on the righthand 
side of (3.8) is at least one degree lower than the operator p (E) on the 
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left-hand side. This can readily be seen by using the definitions (2,3) and 
writing out the highest degree term of the operator 


(3.10) 


o(E) - mp(E) = 




Consequently, the right-hand side of (3.8) can be computed explicitly, that 
is, from known data when advancing the numerical solution from n+k-1 to 
n+k. 


Finally, to apply the combined scheme (3.8) to the PDE (1.1), we split the 
linear differential operator of (3.1), i.e., 

(3.11) f,(u) = (a ——r + c — r^u and f 2 (u) = b - 9 U . 

1 \ 3x 2 3y 2 / 2 8x 3y 

By substituting (3.11) into (3.8), we obtain 

(3.12) f-»At(a^ + c^]p ( E>u" 

* 4t ( a -^2 + C '^') [0(E) ■ “PC E >3“ n + Atb o e (E)u n . 

This formula is implicit for u v „ and u_ T and explicit for u^ r . 

XX yy xy 

Remark ; The scheme (3.12) with the simplest evaluation of the right- 
hand side has a e (E) given by 

(3.13) a e (E) = [a(E) - u>p(E)]u n , 


in which case 

( ' 2 2 2 \ 

a + b t— r — + c -M[c(E) - mp(E)]u n . 

3x 2 8x3y 3y 2 J L 

Unfortunately, the method with cr e (E) defined by (3.13) is only first-order 
accurate. 

For the stability analysis of (3.12), we consider the second-order, two- 
step methods where p(E) and a(E) are defined by (2.7a) and (2.9) and the 
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explicit operator o e (E) is obtained from (2,9) by setting 8=0, 

(3.15) o e (E) = ^ + |^E - (t + j 

The details of the stability analysis are given in Appendix A. Scheme (3.12) 
is found to be unconditionally stable for all values of a,b,c satisfying 
equalities (1.1c, d) if and only if 

(3. 16a, b) 5 < 8 - 1 , 5 > . 

The values of the parameters (8,5) satisfying these inequalities are shown by 
the shaded region of Fig. 2. Inequality (3.16a) is more restrictive than 
(2.16a) for the generating second-order, two-step method to be A^-stable 
(see Fig. 1). Methods that fall in the region between the lines 
(3.17) 5 = 20 - 1 and 5=0-1 

and above 5 = - 1/2 are not unconditionally stable for all values of the 
coefficient b satisfying inequality (l.ld). Note that, with the exception 
of the two-step trapezoidal formula, none of the methods listed in Table 1 
falls in the shaded region of Fig. 2. 

Remark: It is of interest to note that the explicit operator (3.15) can 
also be obtained from the implicit operator (2.9) applied to f n , 

o(E)f n = 8f n+2 + (5 - 20 + |)f n+1 - (5 - 6 + -|)f n , 

by using linear extrapolation, that is, 

f n+2 = 2f n+1 - f n + 0(At 2 ) , 

^ | n 

to approximate f . Although the use of linear extrapolation might sound 
rather disreputable, the application of an explicit LMM does not. 

In the following section we find the rather remarkable result that an 
approximate factorization of the left-hand side of (3.12) into a product of 
one-dimensional operators restores most of the parameter space (0,5) for 
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unconditional stability lost by the explicit treatment of the mixed 
derivative. 

4. An ADI scheme: the p (E) or A formulation. In the procedure sug- 

gested in [20] for designing ADI schemes, the implicit operator to be inverted 
is constructed so that the unknown variable to be determined at each time step 
is p(E)u n . In the absence of a mixed derivative this choice ensures that 
approximate factorization into a product of one-dimensional operators does 
not upset either the temporal accuracy (second-order) or the unconditional 
stability of the scheme. In this section we construct an ADI scheme in the 
p (E) formulation with the mixed derivative treated explicitly with second- 
order accuracy and examine the stability of the resulting algorithm. 

If the spatial derivatives appearing in (3.12) are replaced by appropriate 
difference quotients, then one obtains in general an enormous linear system to 
solve for p(E)u n . This difficulty can be overcome by an approximate factor- 
ization of the left-hand side of (3.12) which reduces the problem to a product 
of one-dimensional spatial operators; that is, 

(4.1) ^1 - wAta - to Ate -^-^p(E)u n 

= At ^a ^2 + c ~ 2 ^[o(E) - o) P (E)]u n + Atb o e ( E) U n . 

Comparison of the left-hand sides of (4.1) and (3.12) shows that they differ 
by the cross-product term 

(4.2) w 2 At 2 a c -^-r p(E)u n . 

dx z 3y 2 

But by expanding p(E)u n in a Taylor series about u 11 and using the consis- 
tency and normalization conditions (2.5a,b), one obtains 
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(4.3) 


a > 1 . 


p(E)u n = At + 0(At“ +1 ) , 


Consequently, the cross-product term (4.2) 


3 2 3 2 n 3 2 3 2 3 n , 

rn^At^a x c — p (E)u = (o^At^a — — c — u + 0(At H ) 


(4.4) 


3x 


3y^ 


3x 2 3y 2 9t 


= 0(At 3 ) 


is a third-order term and the formal accuracy of the scheme (3.12) is not 
upset by the approximate factorization (4.1). 

In practice, the second-order, two-step methods defined by (2.7a) and 
(2.9) are of primary interest and in the remainder of this section we li m it 
our attention to these methods. The explicit operator cr e (E) is given by 
(3.15). In numerical algorithms for partial differential equations it is 
conventional to use n+1 as the most advanced time level; hence, we multiply 
(4.1) by E -1 . The shifted difference operators for the second-order, two- 
step method can be written as 

(4.5a) Au 11 = E _1 p(E)u n = [(1 + 5)E - (1 + 25) + 5E _1 ]u n , 

(4.5b) = [(1 + 5)A - 5V]u n , 


(4.5c) 

= (1 + 5)u n+1 - (1 + 25 )u n + 5u n_1 , 

(4. 6a) 

E 1 o(E)u n = |eE + ^-20+|) 

- ( 6 - 6 + i)e 

(4.6b) 

= |l + 0A + ^ - 0 + 

t)v] u n , 

(4.7a) 

E‘ 1 a e (E)u n = ^5 + |) - (5 + j) I 

1 J" ’ 

(4.7b) 

- [l + (C + i)v]u n . 



where the symbols A and V are classical forward and backward difference 
operators defined by 
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(4.8) 


. n n+1 n 
Au = u - u 


_ n n n-1 
Vu = u - u 


As a notational convenience, we have denoted the operator E ^p(E) by A. 
From (4.5) and (4.6) there follows 

T7~ 1 


(4.9a) 


fa(E) - wp(E)]u n = - w + .2^ - ^ - a) + -0E ^ u n 

= |l + ^ - a) + f)vju n , 


(4.9b) 
where 
(4.10) 

The factored scheme (4.1) becomes 


m = 


1 + 5 ' 


(4.11) (l - wAta r ) f 1 - mAtc •— rjAu 11 

\ 3x 2 /V 8y 2 / 


= At [a 


( a 


3x z 


+ c 


3y 


■)[ + («-» + |)v]u n + 4tb ^ [i + (s + |)v]u” . 


The computational sequence to implement the factored scheme (4.11) as an ADI 
method is not unique, but an obvious choice is 


(4.12a) 

(■- 

wAta JAu* = RHS (4 

3 x 2 / 

(4.12b) 

(- 

0)AtC -^-r-^Au n = Au* , 
3y 2 / 

(4.12c) 


(1 + S)u n+1 = Au 11 + 

where Au* 

is a dummy temporal difference 


n 


n- 1 


Remark: If the first-order explicit method with 0 e (E) defined by (3.13) 

were used rather than a second-order method, the right-hand side of (4.1) 
would be the same as (3.14). In this case, the right-hand side of (4.11) 
would be replaced by 


(4.13) 


At 


( 

\ 3x 2 


, , 3 2 , 

+ b 3^7 + C 


^)[l + ( 5 - 0J + !),]„" 
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where we have used (4.9). Although the resulting scheme is first-order accu- 
rate in time for the mixed derivative, it is unconditionally stable for values 
of (6,5) in the shaded region of Fig. 1. This result was established in [20, 
Appendix A] by constructing the scheme so that the characteristic equation 
(which determines the stability of the method) for the factored partial dif- 
ference equation had the same form as the characteristic equation (2.13) for 
the ordinary differential equation. 

The scheme (4.11) is second-order accurate in the mixed derivative but 
the characteristic equation no longer has the form (2.13). Consequently, the 
stability analysis must be redone in a less elegant manner than in [20] and is 
carried out in Appendix B. The stability analysis of Appendix B is for the 
general two-step ADI scheme developed in the following section. The A for- 
mulation (4.11) is a special case of the general two-step scheme and is found 
to be unconditionally stable for all values of a,b,c satisfying inequalities 
(1.1c, d) if and only if 


(4 . 14a, b) 


0 > 


2(1 + g ) 2 

3 + 4g ’ 



The parameter space (0,5) satisfying these inequalities is shown by the shaded 
region of Fig. 3. The inequality (4.14a) is more stringent than (2.16a), and 
methods that fall in the region between the curves 

(4.15) 5 = 20 - 1 and 0 = - 2(1 + p 2 

3 + 45 

and above 5 = - 1/2 are not unconditionally stable for all values of the 
coefficient b satisfying inequality (l.ld). This includes such popular 
schemes as the trapezoidal formula and Lees method (see, e.g., [2]). However, 
there remains a large class of unconditionally stable methods including, e.g., 
the backward differentiation formula (see Table 1). 
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5. A general two-step ADI scheme. In the ADI formulation described in 


the previous section it is essential that the unknown variable be at least a 
first difference to ensure that the approximate factorization does not degrade 
the second-order accuracy of the scheme; for example, the unknown variable 
E -1 p(E)u n = Au 11 is an approximation to At(8u/9t) (see Eq. (4.3)). The choice 
of unknown variable is not unique and in fact the most general form of a first 
difference using data from the three levels n-1, n, n+1 can be written as 

(5.1) u n+1 - (1 + a)u n + au 11 1 = (A - aV)u n 

= (1 - a)At + (1 + a) — y- + 0(At 3 ) , 

where A and V are the forward and backward operators defined by (4.8) and a 
is an arbitrary real constant. The undivided difference (5.1) is a second 
difference if a = 1 and a first difference otherwise. The most accurate 
first difference is for a = -1. 

5.1. General formulation. A general formulation of two-step ADI schemes 
is obtained if we choose (A - aV)u n as the unknown variable. We first rewrite 
the LMM (3.7) in a convenient form for the construction of an ADI scheme with 
(A - aV)u n as the unknown variable. Multiplying (3.7) by E -1 and inserting 
the operators E -1 p(E), E _1 o(E), and E -1 a e (E) defined by (4.5b), (4.6b), and 
(4.7b), one can rewrite the resulting two-step method as 

(5.2) Au 11 - rnAtAf j 11 = |l + ^ - 6 + f)vj(f? + fj) 

+ rnAtVfz + Y~+T VuT1 

After subtracting aVu 11 - awAtVf] 1 from both sides, there follows 
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(5.3) (A - aV)u n - u>At (A - aV)f? = |l + - 6 + ^ vj ^ 

?2 + aaiAtVf^ + ^ ^ ~ Vu 11 . 


+ wAtVf, 


Using (3.11), one obtains 


(5 


. 4 ) 1 - wA t ^a 


3x^ 


+ c 


JL\\ 

9y 2 /_ 


(A - aV)u 


n 


-i^( a 


+ wAtb 


3x^ 


+ b 


n 


3x3y 


+ c 


9y' 


3x3y 


Vu + awAt (a 


( a 


3x" 


:)[ l + (f - 6 + f) , ] u " 

+ 0 '$) , “ n + (r^r - “) ,un 


The spatially factored form of (5.4) is 


(5.5) 


(i - Mta ^)(l - »Atc 


(A - aV)u n = RHS(5.4) , 


which can 
(5. 6a) 

(5.6b) 


be implemented as 



3 2 

toAta — - )(A 


3x' 


uAtc I (A 
3y 2 - 


> 

* 


- ctV)u 


- aV)u 


n 


RHS (5. 4) , 

(A - aV)u* , 


(5.6c) u n+1 = (A - aV)u n + (1 + a)u n - au 11 * . 

The stability analysis for the general factored scheme (5.5) is given in 
Appendix B. The scheme is found to be unconditionally stable for all values 
of a,b,c satisfying inequalities (lc,d) if and only if 


(5.7a,b,c) 0 > 


2(1 + V) 


2 + 


V 


(1 + a)(l + 25) 
1 + 5 



-1 < a < 1 . 


The parameter space (0,5) satisfying these inequalities is indicated in Fig. 4 
for several values of a. For a given value of a, the stable range is to the 
right of the curve labeled with that value of a and above the curve 
5 = - 1/2. The extent of the (9,5) parameter space for unconditional 
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stability is a monotone increasing function of a in the range [-1,1] with 
the smallest region for a = -1 and the largest for a = +1. Inequality (5.7a) 
is more stringent than (2.16a) as is apparent in Fig. 4. Along the line 
5 = 0, inequality (5.7a) becomes 

(5.8) e > , (5 = o), 

2 + /I + a 

and the smallest allowed value for 0 occurs when a = 1, in which case (5.8) 
becomes 

(5.9) 0 > - 0.586 , (K = 0). 

1 + /I72 V ' 

Along the lower stability boundary 5 = - 1/2, inequality (5.7a) becomes 

(5.10) 9 > | , = 

which is independent of the parameter a. As a consequence of inequali- 
ties (5.9) and (5.10), such popular implicit methods as the trapezoidal for- 
mula (0 = 1/2, 5=0) and Lees method (6 = 1/3, 5 = - 1/2) are not uncondi- 
tionally stable for all values of the coefficient b satisfying 
inequality (l.ld). 

If the parameter a is chosen to be 5/(1 + 5)» scheme (5.5) reduces to 
the A formulation of Sec. 4. This ADI method has the peculiar property that 
the unknown variable depends on the parameter 5> that is, on the particular 
LMM chosen. The parameter space (0,5) for which the A formulation is uncon- 
ditionally stable is shown by the shaded region of Fig. 2 and the extent of 
the region is nearly as large as that for (5.5) with a = 1. 

5.2, Special cases of general formulation. Various constant values of 
a in the range [-1,1] produce useful and interesting algorithms and we con- 
sider several in greater detail. The schemes are named according to the 
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classical difference operator represented by (5.1) for the chosen values of 
a. 


5.2a. The A formulation (a = 0). If a is chosen to be zero, the 
general scheme (5.5) reduces to 


(5.11) 


(1 - wAta - — |(l - oiAtc -^-—)Au n 

V 8y 2 / 

At ( _l£ 

1 + 5 \ a 9 x 2 


+ b 3x3y + C 9y 2 


)[ 1 + (f - 6 + i) 7 ] u " 


2 

+ uAtb 3^7 Vu " + TTJ VuU ’ 

which we call the A formulation [1,19] since Au 11 is the unknown variable. 
The parameter space for unconditional stability is given by the inequali- 
ties (5.7a,b) with a = 0 and is shown graphically in Fig. 4. 


5.2b. The 5 2 formulation (a = 1) . In the general ADI formulation (5.5) 
the unknown variable (A - aV)u n is an approximation to At(3u/3t) if a ^ 1 
(see Eq. (5.1)). A less natural choice for the unknown variable for a first- 
order (temporal) differential equation is the second difference obtained when 
a = 1. In this case, (5.1) becomes 

(5.12) (A - V)u n = <S 2 u n = At 2 - + 0(At 4 ) , 

3t 2 

where the classical second difference operator 5 2 is defined by 

/r i o\ o ? ^ nH - 1 »-> n , n 1 » 

(5. 13) O^u = u - 2u + u 

The possibility of using 6 2 u n as the unknown variable does not arise for 
linear one-step methods (e.g., the trapezoidal formula (2.11)) since these 
methods only involve two time levels. 
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If a 


is set equal to one in the general two-step scheme (5.5), we 


obtain 

(5.14) 


-T*h + b + {t + i}} n - rr 


Vu 


3x z «dy 3v z / 1 \ A / I J- T 5 

A possible advantage of the 6 2 formulation is that the cross-term error 
introduced by the approximate factorization is one order higher than in the 
general (A - aV) formulation with a ^ 1. By comparing the left-hand sides 
of (5.4) and (5.5), we find they differ by 

w 2 At 2 a d ^ c — (a - aV)u n , 


3x z 3y^ 


which, for a = 1, becomes 
(5.15) 


2 At 2 a c 6 2 u n = u> 2 At 4 a c -^- — 
8x 2 3y 2 3x 2 3y 2 3t 2 


+ 0(At 5 ) , 


= 0(At 4 ) , 

where we have used (5.12). The cross-term error for the A = E _1 p(E) formula- 
tion is given by (4.4). The parameter space for which the S 2 formulation is 
unconditionally stable is given by inequalities (5.7a,b) with a = 1 and is 
shown graphically in Fig. 4. 

A distinct disadvantage of the 6 2 formulation occurs when it is applied 
to convective (hyperbolic) model equations as briefly discussed in Section 8. 


5.2c. The 2y£ formulation (a = -1). As a final special case of the 
generalized formulation we choose a = -1, that is, 

(5.16) A - aV = A + V = 2p6 , 
where 2y6 is the classical central difference operator 

(5.17) 2y<5u n = u n+1 - u 11-1 . 
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The scheme (5.5) becomes 


(5.18) ^1 - wAta ~ “Ate _ ^'^2v6u n 


(a — + b r~ ~ + c 


5 - 20 + i)v 

\ 3x 2 3x3 F 3y 2 / 

V 

2/ _ 


n 


+ 2wAtb 


n n , 1 + 2( n n 
Vu + — — : — r 2 - Vu 


3x9y 1+5 

The parameter space for unconditional stability of this formulation is given 
by inequalities (5.7a,b) with a = -1, 


(5.19) 


6 > 1 + 5 , 


5 *-2 ’ 


which is identical to the stability space for the unfactored scheme (3.12) 
(see inequality (3.16) and Fig. 2). 


5.3. Algorithm selection. From the class of unconditionally stable 
methods one can choose a scheme with properties that are desirable with regard 
to computer storage, computational simplicity, and temporal behavior when 
applied to stiff problems and/or problems with nonsmooth data. The choice 
generally requires a compromise. 

Consider, for example, the A formulation (4.11) of Sec. 4, The compu- 
tation of the right-hand side of (4.11) is obviously simplified if we set 

(5.20) 5 - u> + y = 0 . 

Since w = 9/(1 + 5), this equation can be rewritten as 

(5.21) 9 *= <5 + l)(s + |) ; 

it is plotted in Fig. 3. Another variant is obtained by rewriting the right- 
hand side of (4.11) as 
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(5.22) RHS(4.11) = At^a + c ^ - uj + - to + ^E _1 Ju n " 

+ ^ 3S7 [■ + ( e + i}y • 

where we have used (4.9a). The calculation of (5.22) is simplified if we let 

(5.23) K - id + -| = 0 , 


which can be rewritten as 

(5.24) 0 = U + 1)(? + |) , 

it is also plotted in Fig. 3. If in addition, we choose E, = - 1/2, then 
(5.22) becomes simply 

/ 2 2 \ 2 n 

(5.25) RHS (4 . 11) = At (a + c -^Ju 11-1 + Atb . 

\ 3x z 3y 2 / 3x3y 


In this special case, each spatial derivative on the right-hand side of (4.11) 
requires evaluation at only a single time level. The time differencing 
(0 = 1/2, 5 = - 1/2) corresponds to the two-step trapezoidal formula (see 
Table 1). A Q - and A-stable methods along the bottom boundary Z = - 1/2 of 
Fig. 1 are "symmetric" schemes. These methods have the unfortunate property 
that the modulus of at least one root of the characteristic equation (2.13) 
approaches 1 as XAt -* °°. Consequently, these methods can produce slowly 
decaying numerical oscillations when applied to stiff problems. This 
observation illustrates that computational simplicity should not provide the 
sole basis for selecting a time-differencing scheme. 

The computation of the right-hand side of the A formulation (5.11) is 
obviously simplified if we set 

(5.26) S - 0 + j = 0 . 

This special case of the A formulation was given by the authors in [1; see 
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Eq. (4.2)]. The A formulation is particularly attractive for its simplicity 
in programming logic and minimal computing storage requirements. Finally, the 
6 2 u formulation (5.14) has the computational advantage that all spatial 
derivatives on the right-hand side operate on the same function, that is, 

[1 + (£ + l/2)V]u n . For the special case £ = - 1/2, the function is con- 
veniently u 11 . 

Although considerations of computer storage and computational simplicity 
may not be particularly important for the simple model equation (1.1), they 
are of primary concern when one deals with more complicated parabolic equa- 
tions such as the compressible Navier-Stokes equations (see, e.g,, [1,21]). 


5.4. General formulation with no mixed derivative. It is important to 
note that if b = 0, that is, there is no mixed derivative, then inequali- 
ties (5.7) are replaced by 


(5.27) 


£ < 20 - 1 , £ > --| , -1 < a < 1 , 


and the general two-step ADI formula (5.5) is unconditionally stable for the 
same values of (0,£) as for the original second-order, two-step method (see 
inequalities (2.16) and Fig. 1). In the absence of mixed derivatives, the 
natural extension of (5.5) to three spatial dimensions is also unconditionally 
stable for values of (a , 0,£) satisfying inequalities (5.27). 

It is appropriate to mention the relation between the Douglas-Gunn method 
[9, Sec. 3] for multilevel difference schemes and the general two-step ADI 
scheme (5.6) in the absence of a mixed derivative, that is, b = 0. The dif- 
ference (5.1) corresponds to the difference 

(5.28a) 

in the Douglas-Gunn development where 


n+i n+l 
u " u. 

k 
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(5.28b) 


n+1 , n , , n-1 , , , _ , 

u * = <p o u + ^1 U > 4>o + ^ 1 - 1 * 

Hence, on comparing (5.28) and (5.1), one finds that 
(5.29) <J> 0 = 1 + a , 4> 1 = -a . 

Douglas and Gunn give a formal procedure for devising an ADI scheme from a 
fully implicit scheme supplied by the user. For example, consider the second- 
order, two-step method ((3.2) with b = 0) , where p (E) and o(E) are defined 
by (2.7a) and (2.9) and (0,5) satisfy inequalities (2,16). If we apply the 
formulas (3.7) of [9], we obtain an ADI algorithm that can be shown to be 
equivalent to (5.6). The resulting scheme is unconditionally stable for 
-1 < a < 1 since the LMM is A 0 ~stable. Recall that the discussion of this 
paragraph is only for the case of no mixed derivative. 

6. Time-dependent coefficients. If the coefficients a,b,c of the PDE 

(1.1) are functions of time, a difficulty arises when we insert (3.11) into 
(3.8) since 

(6.1) p(E)[a(t n ) + c(t n ) ^-"U[a(t n ) ~^~2 + c(tI1) -^-]p(E)u n . 

L 9x z 9y 2 J L 9x 2 3y 2 J 

This is not an equality because the time dependence of the coefficients cannot 
be neglected when the temporal-difference operator p(E) is applied. This 
problem can be avoided if we begin with the one-leg method (2.21) instead of 
the conventional LMM. 

With Oj and o 2 defined by (3.6), the one-leg method (2,21) is 

(6.2) p(E)u n = Atf (cr(E)u n ,a(E)t n ) + Atf 2 (a e (E)u n ,a e (E)t n ) . 

For the PDE (1.1) we identify f 1 and f 2 as (3.11) and obtain 
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(6.3) 


(E)u 11 = At a(t) — r + c(t) 


3x" 


3 2 *] /T.N-n 

— y a(E)u 

sy 2 J 


+ “[ b < E e> 3^]a e CE>a n 


where t and t e are defined by 


(6. 4a, b) 


t = a(E)t n , t e = a e (E)t 


n 


If (6.3) is modified to the prefactored form by subtracting 


. c 


mA 1 1 a ( t ) — r + c ( t ) 


3x" 


t) p(E)u n 

3y J 


from each side, we obtain 


(6.5) 


[■ 


rTx 3‘ 




1 - o)At a(t) — p + c(t) „ 
3x 2 3y 2 J 


•p(E)u 


[• 


2. 2 

= At | a(t) + c(t) ~~r | [p (E) - mp(E)]u l 
3x 2 3y 2 


■) 


+ Atb (t„) 


Pp(E)u . 


'e' 3x3y “e' 

The prefactored form (6.5) is identical to (3.12) where a, c, and b are 
evaluated at t and t e defined by (6.4). Consequently, the factored scheme 
(4.1) is valid for time-dependent coefficients provided a,c,b are evaluated 
at the appropriate times t and t e . 

For second-order, two-step methods, the shifted-dif ference operators are 
defined by (4.6) and (4.7). For this case 


(6. 6a) 
(6. 6b) 


E -1 t = E _1 a(E)t n = t n + (z + -|^At , 
E“ 1 t £ = E“ 1 a e (E)t n = t n + ^ + yjAt 


and hence the time-dependent coefficients a,b,c are all evaluated at the 
same time which we denote by 
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(6.7) 


t n + (5 + -|^At . 


Therefore, the ADI scheme (4.12) is valid for time-dependent coefficients 
evaluated at t*. Likewise, the same statement applies to the general two- 
step ADI scheme (5.6). 

7. Numerical examples. In this section the ADI methods of the previous 
sections are used to solve the parabolic equation (1.1) for a test problem 
with variable coefficients. The purpose of these numerical experiments is not 
to find the optimum scheme but to demonstrate by numerical example that each 
of the formulations — A, A, and 6 2 — achieves the purported second-order 
accuracy. In addition, we demonstrate the detrimental effect on the accuracy 
if the mixed derivative is treated with a first-order-accurate method or the 
variable coefficients are not evaluated at the proper time level. 

For the example problem, the coefficients are 

(7.1a) a(x,y, t) = x 2 + y 2 ^ (1 + t 2 ) , 

(7.1b) b (x,y , t) = -(x 2 + y 2 ) (1 + t 2 ) , 

(7.1c) c (x, y , t ) = ^x 2 + y 2 ^ (1 + t 2 ) . 

An exact solution is 

(7.2) u(x,y,t) = (x 2 y + xy 2 )exp£-^l + -^-)tj 

Numerical solutions were computed on the unit square (0 < x, y < 1) with the 
initial and boundary values computed from (7.2). For example, the initial 
condition is 

(7.3) u(x,y,0) = x 2 y + xy 2 , 0 < x, y < 1 . 

This model problem is a variant of an example given by McKee and Mitchell [15] 
modified so that the coefficients a,b,c are time-dependent. 
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In the numerical computations of this section, the spatial derivatives 
were approximated by the follox^ing central difference approximations 


(7.4) 


(7.5) 


(7.6) 


9 2 Q 


9x z 

3 2 Q 

9y 2 

5Q 


« 2 .Q 4 


j >k 


w X • 1 

- + 0(Ax 2 ) , 

Ax 2 


9x9y 


j >k 


3 *k 


<5 2 Q- , 
y 3 >k 

Ay 2 


+ 0(Ay 2 ) , 


(pS) x (u6) y 

AxAy 


Q. + 0(Ax 2 ,Ay 2 ) 
1 > k 


(Q, 


- Q,. 


~ Q-s 


+ Q 4 


,) > 


4AxAy y j + l,k+l v j + l,k-l j-l,k+l x j-l,k-l 

where x = jAx and y = kAy. Here 6 and p are classical finite-difference 
operators defined by 

S X«j ‘ Vl/2 - Vl/2 ’ 2 % Q j ‘ q i + !/2 + ^ j — 1 / 2 ’ 

and hence 

5 x Q j = V> - 2Q i + Q i-1 ’ 

2(p6) x Q j = Q j+1 " Q j-l’ etC * 

Consider, for example, the A formulation (4.12). With the spatial 
derivatives replaced by the central-difference quotients (7.4)-(7.6), the 
x- and y-operators on the left-hand side of (4.12a,b) each requires the solu- 
tion of a tridiagonal system. There is a well-known and highly efficient 
solution algorithm for tridiagonal systems (see, e.g., [11, p. 551). The 
solution of the x-operator (4.12a) (along each y-constant line) requires 
the dummy temporal difference Au* along the left and right boundaries. In 
problems considered in this section, we assume that u(t) is given on the 
boundaries, and consequently Au* can be determined by an explicit calcula- 
tion using (4.12b) applied along both the left- and right-hand boundaries. 
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This is the initial calculation made when advancing the solution from n to 
n+1. Application of the general two-step ADI scheme (5.6) requires an analo- 
gous computation of a dummy temporal difference along the left and right 
boundaries. 

Since the algorithms considered in this paper are, in general, two-step 
(temporal) schemes, a solution at t = At is needed together with the initial 
condition to start the computation. For the numerical examples computed 
herein, the exact solution (7.2) at t = At was used to provide the addi- 
tional level of data. In practice, one can use (4.12) as a one-step method on 
the first time step. This is accomplished by replacing the right-hand side of 
(4.11) by (4.13) and choosing 0=1/2, £ = 0. 

The numerical differentiation formulas (7.4) and (7.5) are exact (i.e., 
the truncation error is zero) for a polynomial of degree not exceeding three 
and (7.6) is exact for a polynomial of degree not exceeding two. Since the 
exact solution (7.2) is a quadratic polynomial of degree two in each spatial 
variable, the numerical solution of an unfactored algorithm would have the 
peculiar property that there would be no spatial discretization error. Con- 
sequently, the error in a numerical solution for the example problem (7.1) 
consists of the temporal discretization error and the cross-product error term 
from the approximate factorization (see, e.g. , Eq. (4.4)), and, of course, 
roundoff error. 

For each numerical experiment, we compute the L 2 norm of the error 

which is defined as follows. At a given time, t = t n = nAt, the error e. , 

J » K 

at each grid point is defined by 

(7.7) e^ k = u”^ - u(jAx, kAy, t n ) , 
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where u° , is the numerical solution and u(jAx,kAy ,t n ) is the analytical 
J > K 

solution. The Euclidean or L 2 norm of the error is defined by 



where J and K are the total number of grid points in the x- and y-directions. 

The second-order backward differentiation method (0=1, £ = 1/2) (see 
Table 1) was chosen as the generating LMM for the first computational experi- 
ment. The L 2 errors for the A, A, and 6 2 formulations (algorithms (4,12), 
(5.11), and (5.14)) are shown in Table 2. Each computation was carried out to 
a given time (t = 1.0) with a fixed ratio of At/Ax = 1.0. The computations 
were repeated with successively smaller values of At so that the L 2 rate 
could be computed. (The L 2 rate is the slope of a log-log graph of the L 2 
error vs. At. For a second-order method without round-off error, the L 2 
rate should approach two as At 0.) The results show the second-order accu- 
racy of the methods. Since the same time-differencing method was used for 
each computation, the differences in the L 2 error (for a given At) result 
from the cross-product error of the approximate factorization. 

The next numerical experiment demonstrates the detrimental effect on the 
accuracy if the mixed derivative is treated with first-order accuracy. The 
errors listed in Table 3 were computed using the A formulation with the back- 
ward differentiation method as the generating LMM. For reference, the results 
listed under column (1) are repeated from Table 2. The L 2 errors and rates 
tabulated under column (2) were obtained using (4.12) but with the right-hand 
side of (4.11) replaced by (4.13), that is, a first-order temporal treatment 
for the mixed derivative. The degradation in accuracy is obvious. 
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Column (3) of Table 3 shows the deterioration in accuracy when the time- 
dependent coefficients a,b,c are not evaluated at the proper time level. 

The coefficients should be evaluated at t* defined by (6.7) and hence for 
the backward differentiation method (5 = 1/2) t* = t n + At. In obtaining the 
L 2 errors listed in column (3), the coefficients were evaluated at t rather 
than t* and the loss of accuracy is apparent. 

It is important to note for 5=0 that the operator A defined by (4.5) 
in the A formulation becomes 

, n . n 
Au = Au 

Consequently, the A and A formulations are identical if 5 = 0. (Recall 
that the A algorithm is given by (5.6) with a = 0.) An advantage of the A 
formulation is that u 11 1 is not needed to compute u n+1 in the final step 
(5.6c); hence, the A form generally requires the least amount of storage. 

On the other hand, the A formulation for 5^0 has a significantly 
reduced parameter space (0,5) for unconditional stability when applied to 
hyperbolic problems (see Sec. 8 and [20, Sec. 9]). Consequently, because the 
A and A formulations are identical for 5=0, this subclass of schemes has 
the simplicity of the A form and the robustness of the A form. Table 4 
compares the L 2 error and rate for several schemes with 5 = 0 * For a 
fixed value of 5 in the region of A 0 -stability (see Fig. 1), the error con- 
stant is a monotone increasing function of 8. This is verified by comparing 
the L 2 errors for a given value of At in Table 4. 

According to the stability analysis of Appendix B, methods that fall in 
the region between the curves (4.15) are not unconditionally stable in the A 
formulation (4.11) for all values of the coefficient b satisfying inequality 
(l.ld). The last numerical experiment verifies this result for the LMM 
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methods listed in Table 1. The parameters used in the computation are listed 
in the caption of Table 5. The loss of unconditional stability for the trape- 
zoidal formula and Lees method is apparent from the large error for these two 
methods listed in Table 5. 

8. Concluding remarks. In this paper we combine the class of all 
A Q -stable second-order linear two-step methods and the method of approximate 
factorization to construct unconditionally stable ADI methods for parabolic 
equations (in two space dimensions) with a mixed derivative. For computational 
simplicity the mixed derivative is treated explicity. 

In the application of the approximate factorization method we consider 
several different solution variables. In Sec. 4 we follow [20] and select 
p(E)u n as the unknown variable. This choice provides a natural framework for 
constructing unconditionally stable ADI methods for parabolic PDEs by combin- 
ing A Q -stable LMMs with approximate factorization. The choice of the unknown 
variable is not unique and for completeness we derive a general two-step ADI 
scheme with (A - aV)u n as the unknown variable in Sec. 5. The general formu- 
lation contains a parameter a in addition to the parameters (0,£) of the 
second-order, two-step method. The parameter space (a, 0,0 for unconditional 
stability is determined in Appendix B. Several general observations can be 
made regarding the stability of these schemes: (1) For a given value of a 

in the range [-1,1], the parameter space (0,C) for unconditional stability is 
reduced from that of the unfactored implicit algorithm (3.2) (compare Fig. 1 
with Figs. 3 and 4), but is increased from that of the unfactored implicit- 
explicit (i.e., explicit treatment of mixed derivative) algorithm (3.12) 
(compare Fig. 2 with Figs. 3 and 4). (2) For any allowed value of a, the 

reduced parameter space excludes the familiar (one-step) trapezoidal formula 
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and the Lees type scheme (see Table 1). (3) The 6 2 u formulation retains the 

largest parameter space for unconditional stability. (4) The p (E) or A 
formulation has the peculiar property that ct = 5/(1 + 5), that is, a depends 
on the particular LMM chosen. The extent of the parameter space for uncondi- 
tional stability (Fig. 3) is nearly as large as for the 6 2 formulation 
(Fig. 4 with a = 1). (5) Although not considered in this paper, it can be 

shown that if the general (A - aV) formulation is applied to the model convec- 
tion equation 


( 8 . 1 ) 


3u 


3u 


3u 


3t C 1 3x C 2 3y » 


then only the p(E) formulation (i.e. , a = 5/(1 + 5) retains the same param- 
eter space (shaded region of Fig. 1) for unconditional stability as the gen- 
erating A-stable LMM. In fact, the 6 2 formulation has no parameter values 
(0,£) for which the scheme is unconditionally stable. 

The emphasis of this paper is on the construction of unconditionally 
stable second-order accurate ADI methods for the model parabolic equation (1.1). 
By following the approach outlined herein (i.e., the use of A Q -stable LMMs 
in conjunction with the method of approximate factorization) one can easily 
construct algorithms for multidimensional nonlinear parabolic systems. For 
some auspicious reason, the parameter space (0,5) for which the class of 
second-order two-step methods is A Q -stable happens to coincide with the 
parameter space for which this class of methods is A-stable. Consequently, 
one can use the class of time-differencing schemes of this paper to design 
second-order ADI algorithms for mixed hyperbolic-parabolic systems of nonlinear 
equations. A noniterative algorithm in the A-form for nonlinear systems was 
considered in [1] and a general development for the p (E) formulation is in a 
companion paper [21]. 
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Appendix A. Stability analysis of combined LMMs. In this appendix we 
examine the stability of the combined LMM (3.7) applied to the model split 
ODE: 

(A.la,b) -^r = X x u + X 2 u ; Xj < 0, X! + X 2 ^ 0 > 

where X2 and X 2 are real constants. In addition, we investigate the stabil- 
ity of the unfactored scheme (3.12) for the PDE (1.1). The analysis is for 
the class of all second-order, two-step methods. 

Consider the combined LMM (3.7) where p (E) , a(E), and o e (E) are defined 
by (2.7a), (2.9), and (3.15). If we apply this scheme to the model equa- 
tion (A. 1) with fj = Xju and f 2 = X 2 u, we obtain a difference equation 
whose characteristic equation is 

(A. 2) a 2^ 2 + a l^ + a o = 0 > 

where 

(A. 3a) a 2 = (1 + O - 6XjAt , 

(A. 3b) aj = -(1 + 20 - U - 26 + J-^XjAt - (z + |)x 2 At , 

(A. 3c) a Q = ? + - 0 + IJXjAt + (? + f) A 2 At * 

Equation (A. 2) is a von Neumann polynomial [16], that is, |c| <1, if and only 

if 

(A. 4a) a 0 < a 2 , 

and 

(A. 4b) -(a 2 + a Q ) < aj < a 2 + a Q , 

where without loss of generality a 2 is assumed to be positive. The inequal- 
ities (A. lb) and (A. 4) lead to the following conditions for the stability of 
the combined two-step scheme: 
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(A. 5a, b) 5 > > 0 < (1 + 2?) + (1 - 20 + 5)A 1 At + (1 + 5)A 2 At . 

In particular, the conditions for A Q -stability are 
(A. 6a, b) 5 ^ 2 » ^2 > ^ ^1 * 

Note that for the special case X 2 = 0> (A. 6b) becomes £< 2B - 1 and condi- 
tions (A. 6) reduce to (2.16). 

For the stability analysis of the unfactored scheme (3.12) for the PDE 
(1.1) we need only consider the stability of the linear two-step scheme (3.7) 
applied to the ODE (3.4) for the Fourier coefficient. In this appendix we 
consider only the spatially continuous solution; however, the results are 
applicable to the spatially discrete case (see last paragraph of Appendix B) . 
The conditions on the parameters (0,5) for the unconditional stability of 
(3.12) can be derived from the A Q -stability requirements (A. 6) and the rela- 
tions 

(A. 7a, b) Aj = -(aKj 2 + ck 2 2 ) , A 2 = -biCjKg > 

obtained by comparing (A.l) and (3.4b). There follows 

(A. 8a, b) 5 > -j , -bK i K 2 - ^ (a<j 2 + c< 2 2 ) . 

The inequalities (A. 8a, b) together with (1.1c, d) imply 
(A. 9a, b) 5 > -j , 5 < 0 - 1 . 

Inequality (A. 9b) is more restrictive than the inequality (2.16a) for the 
generating two-step method (2.6) to be A 0 -stable. Consequently, we have the 
result that the second-order explicit treatment of the mixed derivative 
reduces the parameter space (0,5) for which the unfactored scheme (3.12) is 
unconditionally stable (see Fig. 2). 
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Appendix B. Stability analysis for two-step ADI schemes. In this appen- 


dix we perform a linear stability analysis for the factored scheme (5.5). We 

assume that u n is spatially continuous and seek a solution of the form 

. n n i(K.x + K 0 y) 

(B. 1) u = v e 1 

where v n is the Fourier coefficient and Kj,^ are the Fourier variables. 
Prior to an actual numerical computation, the spatial derivatives are replaced 
by appropriate difference quotients; however, as indicated at the end of this 
appendix, the stability proof for the spatially discrete case requires only a 
minor modification of the following stability proof. 

By substituting (B.l) into (5.5), we find that the Fourier coefficient 
satisfies 


(B. 2) (1 + A)(l + C) (A - aV)v n 


where we have defined 


- (A + B + 


c)[i + (e - e + y)vjv n 

n - a (A + C)Vv n + ~ «)vv n > 


(B. 3) 


A = wAtaKj 2 , B = 0)AtbK 1 ic 2 , C = wAtcK 2 2 


and u - 0/(1 + 5). The amplification factor is defined by 
(B. 4) v n+1 = Cv 11 , 

and consequently it follows from (B.2) that C satisfies the quadratic 
equation 


< B - 5 ) a 2 ? 2 + a 2 c + a Q = 0 , 

where 

(B. 6a) a 2 = (1 + A) (1 + C) , 

aj = - (1 + a) (1 + A) (1 + C) + i U - 9 + I) (A + B + C) 

(B. 6b) 

+ B + a (A + C) - [yTJ ~ a ) ’ 
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(B.6c) 


a Q = o(l + A)(l + C) - | (s - 6 + |)(A + B + C) 

- B - a (A + C) + (j-f-jr - a) . 

In the one-dimensional case (b = c = 0 in Eq. (1.1) and B = C = 0 in 
(B.5)), the roots of the quadratic (B.5) have modulus bounded by unity for 
those values of (0,0 shown in the shaded region of Fig. 1, that is, 

(B. 7a, b) 5 < 26 - 1 , 5 > - j • 

This one-dimensional result follows from the analysis [20]. Note that a 
only enters as a parameter in the two-dimensional factored algorithm (5.5). 
(Recall that (5.2) and (5.3) are actually identical.) One can easily verify 
that the coefficients (B.6) do not depend on a if B = C = 0. We must 
determine if there are additional restrictions on the parameters (8,5) for the 
unconditional stability of the factored scheme (5.5) for arbitrary values of 
a,b,c subject only to the parabolicity conditions 
(B.8a,b) a > 0 , b 2 < 4ac 

of the partial differential equation (1.1). Since the one-dimensional problem 
(b = c = 0) is a special case of the two-dimensional problem, we need not con- 
sider values of (0,5) outside the domain (B.7). Hence to > 0, A and C as 
defined by (B.3) are positive, and 

(B.9) A + B + C =4oAt(aKj 2 + bKji^ + ck 2 2 ) > 0 » 

since the positive definiteness of this quadratic form was the condition which 
led originally to (B.8). 

The coefficients (B.6) of the quadratic (B.5) are real and consequently 
the roots 5 satisfy |?| < 1 if and only if the inequalities (A. 4) of 
Appendix A are satisfied. If we insert a 0 and a 2 as given by (B.6) into 
(A. 4a) , there follows 
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(B. 10) 


0 < + (1 - a)AC + | (? + f)(A + B + C) , 


which is satisfied for all allowable A,B,C if and only if 
(B. 11) a < 1 . 

(Recall that the parameters 0 and E, are required to satisfy inequalities 
(B. 7).) Likewise the left inequality of (A. 4b) is satisfied. If the coeffi- 
cients (B.6) are inserted into the right inequality of (A. 4b) one obtains 

(B. 12 ) 0 < (1 + a)AC - - B + -- + (z - -Va + C) . 

w 1 + £ \ 0)/ 

The determination of necessary and sufficient conditions for this inequality 
is simplified if it is rewritten as 

(B. 13) 0 < (1 + a)ac(k 1 k 2 ) 2 - ^ bkxk 2 + * + ~ 2 ^ - + ^2 - ^(ak^ + ck 2 2 ) 

where we have used the definitions (B.3) and defined: 

(B.l4a,b) k^ = /coAt , k 2 = /uAt k 2 . 

It can be shown that necessary and sufficient conditions for the polynomial P 
in kj,k 2 defined by 

(B.15) P = e J (k 1 k 2 )^ + e 2 kjk 2 + e 3 + e^k^ 2 + e^k 2 2 

to be positive semidefinite (i.e. , P > 0 for all real kj,k 2 ) are 
(B. 16a) e 2 , e 3 , e 4 , e 5 > 0 , 

(B. 16b) |e 2 | < 2/e^eJ + 2/e 4 e 5 , 

Comparison of (B.15) and (B.13) leads to the conditions 

(B. 17a, b) (1 + a)ac > 0 , * 0 , 

(B. 17c , d) ( 2 “ u) a - 0 » ( 2 " w) C - 0 » 

(B.We) |i| < if * c)3c Vrf + 2|/(2 - - • 
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Inequalities (B. 17a,b,c,d) are satisfied by virtue of inequalities (B.7), 
(B.8a), and (B.ll) plus the constraint: 


(B. 18) 


-1 < a . 


Inequality (B.17e) can be rewritten as 


(B. 19) 



2/ac 


V 


(1 + a)(l + 2K) 
1 + K 



which is satisfied for all allowable a,b,c (see inequalities (B.8)) if and 


only if 
(B. 20) 


or 


1 + [(1 +~a) (1 + 2Q 

a) - Y 1 + C + 



(B. 21) 


9 > 


2(1 + K) 


2 + 




(1 + 2Q 
+ l 


Hence the final inequalities which must be satisfied are (B.7b), (B.ll), (B.18), 


and (B. 21) . 


In the above stability analysis we assumed that the spatial derivatives 
were continuous. Since in practice the spatial derivatives are replaced by 
discrete difference quotients, it remains to consider the spatially discrete 
case. If, for example, the spatial derivatives in (5.5) are replaced by the 
second-order difference quotients (7. 4)- (7. 6), then the stability analysis 
proceeds as above with the exception that the exponential in (B.l) is replaced 
by 


(B. 22) 


u n = v n e i(K 1 jAx+K 2 kAy) 
j jk 


where x = jAx, y = kAy. If we make the following correspondence 

2 sin(0 /2) 2 sin(0 2 /2) 

(B. 23a, b) *i t- , < 2 - r~ , 
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(B. 24) 


where 


® 1 

B B cos — cos — 


0! = k ^ Ax , 0 2 = ic^Ay 


between the parameters for the discrete and continuous case, then the amplifi- 
cation factor for the discrete case satisfies the same quadratic (B.5) with 
coefficients (B.6). Since the stability region defined by inequalities (B.7b), 
(B.ll), (B.18), and (B. 21) is valid for arbitrary values of Kj and k 2 and 

01 02 ii 

(B.25) B cos — cos — < | B | , 

we obtain the same stability range for the discrete case. If one uses a non- 
centered approximation for the mixed derivative such as 
9 2 Q 


3x3y 


(B. 26) 


1 

'Vy + W Q j 

2AxAy 

1 

(Q j+l,k " 2Q j,k 

2AxAy 

+ Vo + V-i - Q . 


- Q,. 


where 


A Q, = Q., - Q. , V Q, = Q. - Q, , . etc., 

* J 3+1 J * 3 J J~1 


the only modification necessary in the stability analysis is replacement of 
(B. 24) by 


(B. 27) 


B *■ B cos 
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TABLE 1 


Partial list of second-order two-step methods. 


e 

5 

4> 

Method 

Symbol in 
Fig. 1 

1 

2 

0 

0 

One-step trapezoidal formula 

■ 

1 


0 

Backward differentiation 

• 

1 

3 


1 

"3 

Lees type [14] 

T 

3 

4 

0 

g 

Adams type [17] 

♦ 

1 

2 

1 

~ 2 

1 

Two-step trapezoidal formula 

A 
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TABLE 2 


L 2 error of the A, A, and 6 2 formulations at t = 1.0. 


At = 
Ax = Ay 

Number of 
time steps 

At /Ax 2 

A formulation 

A formulation 

6 2 formulation 

L 2 error 

L 2 rate 

L 2 error 

L 2 rate 

L 2 error 

L 2 rate 

0.2 

5 

5 

0. 785xl0~ 2 

2.02 

2.01 

2.01 

0.107X10 -1 

2.01 

2.03 

2.03 

0.134x10-2 

1.57 

1.72 

1.89 

0.1 

10 

10 

0.193X10" 2 

0.266X10 -2 

0.452xl0" 3 

0.05 

20 

20 

0.479xl0 -3 

0.649xl0 -3 


0.025 

40 

40 

0.119xl0“ 3 

0.160xl0“ 3 

0.372xl0" 4 




































TABLE 3 


Numerical experiments illustrating (1) second-order A formulation, (2) degradation in 
accuracy when mixed derivative is computed with a first-order method, (3) deterioration 


in accuracy when time-dependent coefficients are not evaluated at proper time level. 


At = 

Number of 

At/Ax 2 

(1) 

(2) 

(3) 

Ax = Ay 

time steps 

L 2 error 

L 2 rate 

L 2 error 

L 2 rate 

L 2 error 

L 2 rate 

0.2 

5 

5 

0.758xl0 -2 

2.02 

0.585xl0" 2 

0.59 

0.907xl0- 2 

1.72 

0.1 

10 

10 

0.193X10" 2 

0.389xl0 -2 

0.275X10" 2 

2.01 

0.85 

1.63 

0.05 

20 

20 

0.479X10" 3 

0.216xl0 -2 

0.888X10 -3 

2.01 

0.94 

1.47 

0.025 

40 

40 

0.119X10 -3 

0.113xl0 -2 

0.320X10 -3 















TABLE 4 


L 2 error of the A formulation for £ = 0 and several values of 0 at t = 1.0. 


At = 
Ax = Ay 

Number of 
time steps 

At /Ax 2 

K = 0, 0 = 2/3 


5 = 0, 0 « 3/2 

L 2 error 

L 2 rate 

L 2 error 

L 2 rate 

L 2 error 

L 2 rate 

0.2 

5 

5 

0.705*10“ 2 

2.06 

2.03 

2.02 

0.871xl0“ 2 

2.07 

2.03 

2.02 

0.262X10 -1 

1.85 

2.00 

2.01 

0.1 

10 

10 

0.169xl0 -2 

0.208xl0" 2 

0.726xl0 -2 

0.05 

20 

20 

0.415xl0 -3 

0.511X10 -3 

0.181xl0 -2 

0.025 

40 

40 

0.102X10 -3 

0.126xl0 -3 

0.448xl0 -3 




































TABLE 5 


L 2 error of A formulation (4,11) 
at t = 1.0. Parameters are 
Ax = Ay = 0.025, At = 0.005, number 
of time steps = 200, 


| a . . | At/Ax 2 - 23. 
1 j , k 1 max 


Method 

L 2 error 

One-step trapezoidal 

0.246xl0 2 

Backward differentiation 

0.479*10 -5 

Lees type 

0.405xl0 19 

Adams type 

0.505xl0“ 5 

Two-step trapezoidal 

0.772xl0 -5 






FIGURE CAPTIONS 


FIG. 1. A - and A-stable domain of the parameters (0,5) for the class of all 
second-order two-step methods. Symbols denote methods listed in Table 1. 

FIG. 2. Unconditionally stable domain of the parameters (6,5) for the unfac- 
tored scheme (3.12) with p (E) , a (E) , and cr e (E) defined by (2.7a), (2.9), 
and (3.15). 

FIG. 3. Unconditionally stable domain of the parameters (0,5) for the fac- 
tored A formulation (4.11). 

FIG. 4. Unconditionally stable domain of the parameters (0,5) for the fac- 
tored general formulation (5.5) for several values of a. 
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-1 


% = - 1/2 
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Figure 1.- 
second- 


A and A stable domain of the parameters (9,£) for the class of all 
order two step methods. Symbols denote methods listed in Table 1. 


52 





Figure 2.- Unconditionally stable domain of the parameters (0,5) for the 
unfactored scheme (3.12) with p(E), o(E), and a e (E) defined by (2.7a), 
(2.9), and (3.15). 





Figure 3.- Unconditionally stable domain of the parameters (0,?) for the fac 

tored A formulation (4.11). 
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